Dicom loading and visualization¶

Load DICOM files¶

Libraries¶

In [ ]:
import os
import pydicom
import numpy as np
import matplotlib.pyplot as plt
from typing import List
import cv2
import matplotlib.colors
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation

Aqui alojamos las funciones¶

In [ ]:
def load_dicom_files(directory: str, file_extension: str = '.dcm'):
    """Load DICOM files from a specified directory that match the given file extension."""
    return [pydicom.dcmread(os.path.join(directory, f)) for f in os.listdir(directory) if f.endswith(file_extension)]

def sort_dicom_files_by_position(dicom_files):
    """Ordena archivos DICOM según la posición del paciente en la imagen."""
    return sorted(dicom_files, key=lambda x: float(x.ImagePositionPatient[2]))

def get_pixel_spacing(dicom):
    return dicom.PixelSpacing, dicom.SliceThickness

def sort_dicom_files(dicom_files):
    """Sort DICOM files by available headers, only including those with AcquisitionNumber equal to 2."""
    def get_sort_key(file):
        # Define sorting key: ImagePositionPatient (converted to a sortable form)
        position_patient = getattr(file, 'ImagePositionPatient', None)
        if position_patient:
            position_patient = float(position_patient[2])
        return position_patient if position_patient else 0

    # Filter to include only files with AcquisitionNumber equal to 2 --> With 3D Slicer we have seen is the correct one
    filtered_files = [file for file in dicom_files if getattr(file, 'AcquisitionNumber', None) == '2']
    
    '''
    Primero ordena por AcquisitionNumber y, dentro de cada número de adquisición, 
    ordena por la posición del paciente a lo largo del eje Z (ImagePositionPatient[2]).'''
    # Sort the filtered files by ImagePositionPatient along the Z-axis
    return sorted(filtered_files, key=get_sort_key)


def check_acquisition_uniformity(dicom_files: List[pydicom.Dataset]):
    """Check if all DICOM files have the same 'Acquisition Number'."""
    acquisition_numbers = {img.AcquisitionNumber for img in dicom_files}
    if len(acquisition_numbers) > 1:
        print("Multiple acquisition numbers found. Ensure all CT images are from a single acquisition.")
    else:
        print("Only a single acquisition")
    
def get_pixel_arrays(dicom_files: List[pydicom.Dataset]) -> np.ndarray:
    """Extract pixel arrays from sorted DICOM files and stack them into a single numpy array."""
    return np.stack([file.pixel_array for file in dicom_files])

def apply_cmap(img, cmap_name='bone'):
    """Apply a colormap to an image."""
    cmap = matplotlib.cm.get_cmap(cmap_name)
    normed_data = (img - np.min(img)) / (np.max(img) - np.min(img))
    mapped_data = cmap(normed_data)
    return (mapped_data[:, :, :3] * 255).astype(np.uint8)

def apply_alpha_fusion(dcm_img, seg_img, alpha=0.5, dcm_cmap='bone', seg_cmap='Set1'):
    """Apply alpha blending to two images."""
    dcm_img_colored = apply_cmap(dcm_img, cmap_name=dcm_cmap)
    seg_img_colored = apply_cmap(seg_img, cmap_name=seg_cmap)
    return cv2.addWeighted(dcm_img_colored, 1-alpha, seg_img_colored, alpha, 0)

def display_images(ct_image, seg_image):
    """Display CT and Segmentation images side-by-side."""
    plt.figure(figsize=(12, 6))

    plt.subplot(1, 2, 1)
    plt.imshow(ct_image, cmap='gray', aspect='auto')
    plt.title('CT Slice')

    plt.subplot(1, 2, 2)
    plt.imshow(seg_image, cmap='gray', aspect='auto')
    plt.title('Segmentation Slice')
    plt.show()

def apply_window(image: np.ndarray, window_center: int, window_width: int) -> np.ndarray:
    """Apply window leveling to an image based on the given center and width."""
    lower_bound = window_center - 0.5 * window_width
    upper_bound = window_center + 0.5 * window_width
    return np.clip((image - lower_bound) / window_width, 0, 1)
In [ ]:
# Data paths
# Windows
#ct_directory = "C:/Users/pedromarti/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/4.000000-Recon 2 3 PHASE LIVER ABD-87008"
#seg_directory = "C:/Users/pedromarti/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/300.000000-Segmentation-39839/1-1.dcm"

# Windows feina
ct_directory = "HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/4.000000-Recon 2 LIVER 3 PHASE AP-28011/"
seg_directory = "HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/300.000000-Segmentation-39839/1-1.dcm"
# Mac
#ct_directory = "HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/4.000000-Recon 2 LIVER 3 PHASE AP-28011/"
#seg_directory = "HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/300.000000-Segmentation-39839/1-1.dcm"

Cargamos las imagenes y las ordenamos¶

Cargamos y ordenamos las imagenes CT. Como ya hemos visto en el 3D Slicer, cogemos la acquisition number 2.

In [ ]:
# Load and process DICOM CT images
ct_files = load_dicom_files(ct_directory)
sorted_ct_files = sort_dicom_files(ct_files)

Cargamos la segmentation image

In [ ]:
# Load and process segmentation image
dicom_segmentation = pydicom.dcmread(seg_directory)
segmentation_image = dicom_segmentation.pixel_array

Simple visualización para saber que todo ha ido bien

In [ ]:
# Visualize a slice of CT and Segmentation images
ct_images = get_pixel_arrays(sorted_ct_files)
display_images(ct_images[8], segmentation_image[8])
No description has been provided for this image

Como sabemos que tiene 4 regiones en la segmentación y 300 registros, dividimos 300/4 y así sabemos que cada región tiene 75 registros.

Estudio de encabezados¶

In [ ]:
# Print details and dimensions
print("Details of CT images sorted by Image Position Patient:")
for img in sorted_ct_files:
    print(f"Acquisition Number: {img.AcquisitionNumber}, Image Position Patient: {img.ImagePositionPatient}")

print("Dimensiones del arreglo de segmentación:", segmentation_image.shape)
print("Dimensiones del arreglo de imágenes CT:", ct_images.shape)
Details of CT images sorted by Image Position Patient:
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -240.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -237.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -235.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -232.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -230.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -227.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -225.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -222.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -220.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -217.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -215.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -212.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -210.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -207.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -205.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -202.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -200.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -197.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -195.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -192.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -190.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -187.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -185.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -182.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -180.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -177.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -175.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -172.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -170.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -167.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -165.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -162.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -160.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -157.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -155.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -152.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -150.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -147.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -145.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -142.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -140.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -137.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -135.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -132.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -130.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -127.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -125.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -122.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -120.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -117.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -115.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -112.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -110.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -107.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -105.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -102.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -100.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -97.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -95.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -92.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -90.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -87.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -85.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -82.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -80.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -77.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -75.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -72.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -70.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -67.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -65.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -62.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -60.000000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -57.500000]
Acquisition Number: 2, Image Position Patient: [-192.500000, -180.000000, -55.000000]
Dimensiones del arreglo de segmentación: (300, 512, 512)
Dimensiones del arreglo de imágenes CT: (75, 512, 512)

Ahora ya tenemos la acquisition que mejor queda con nuestra segmentación.

In [ ]:
check_acquisition_uniformity(sorted_ct_files)
Only a single acquisition

Aplicamos windowing a las imagenes para mejorarlas¶

In [ ]:
# Windowing parameters
window_center, window_width = 40, 400
windowed_ct_images = np.array([apply_window(img, window_center, window_width) for img  in ct_images])
In [ ]:
# Tomamos el primer archivo DICOM para obtener los metadatos
first_dicom_file = sorted_ct_files[0]  

# Obtener información del voxel
pixel_spacing = first_dicom_file.PixelSpacing  # Lista: [tamaño del pixel en X, tamaño del pixel en Y]
slice_thickness = first_dicom_file.SliceThickness  # Espesor de la rebanada en direccion Z

# Calcular las proporciones de aspecto
aspect_ratio_coronal = slice_thickness / pixel_spacing[0]
aspect_ratio_sagittal = slice_thickness / pixel_spacing[1]
In [ ]:
import matplotlib.pyplot as plt

def visualize_planes(ct_volume: np.ndarray):
    """Visualize axial, coronal, and sagittal planes of a CT volume with corrected aspect ratios."""
    # Determinar los planos
    axial_plane = ct_volume[ct_volume.shape[0] // 2, :, :]
    coronal_plane = ct_volume[:, :, ct_volume.shape[2] // 2]
    sagittal_plane = ct_volume[:, ct_volume.shape[1] // 2, :]

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # Visualizar el plano Axial
    axes[0].imshow(axial_plane, cmap='gray')
    axes[0].set_title('Axial Plane')
    axes[0].axis('off')

    # Visualizar el plano Coronal con la proporción de aspecto corregida
    axes[1].imshow(coronal_plane, cmap='gray', aspect=aspect_ratio_coronal)
    axes[1].set_title('Coronal Plane')
    axes[1].axis('off')

    # Visualizar el plano Sagital con la proporción de aspecto corregida
    axes[2].imshow(sagittal_plane, cmap='gray', aspect=aspect_ratio_sagittal)
    axes[2].set_title('Sagittal Plane')
    axes[2].axis('off')

    plt.show()

visualize_planes(windowed_ct_images)
No description has been provided for this image

Aplying alpha-fusion¶

In [ ]:
# Cargamos los registros a su region
img_liver = segmentation_image[1:75,:,:]
img_mass = segmentation_image[76:150,:,:]
img_portal_vein = segmentation_image[151:225,:,:]
img_abdominal_aorta = segmentation_image[226:300,:,:]

img_segmented = (img_liver*1) + (img_mass*2) + (img_portal_vein*3) + (img_abdominal_aorta*4)

Visualising the four areas of interest

In [ ]:
import matplotlib.colors as mcolors

# Colormap para las segmentaciones
cmap = mcolors.ListedColormap(['none', 'red', 'green', 'blue', 'yellow'])  # 'none' es para transparencia
bounds = [0, 1, 2, 3, 4, 5]
norm = mcolors.BoundaryNorm(bounds, cmap.N)

def plot_ct_with_segmentations(ct_image, segmentation):
    plt.figure(figsize=(12, 6))

    # Mostrar imagen CT
    plt.imshow(ct_image, cmap='gray')

    # Superponer la segmentación
    plt.imshow(segmentation, cmap=cmap, norm=norm, alpha=0.5)  # Alpha para transparencia

    plt.title('CT Image with Segmented Areas')
    plt.axis('off')
    plt.colorbar(ticks=[1, 2, 3, 4], label='Segments')
    plt.show()

slice_index = 48 # Índice de la slice a visualizar
plot_ct_with_segmentations(windowed_ct_images[slice_index, :, :], img_segmented[slice_index, :, :])
No description has been provided for this image

Prueba: Superposición de las CT sobre la segmentada¶

In [ ]:
# Visualize images and segmentations
for i in range(10):  # Example for the first 10 slices
    plot_ct_with_segmentations(windowed_ct_images[i+5, :, :], img_segmented[i+5, :, :])
print("Segmentation array loaded with shape:", segmentation_image.shape)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
Segmentation array loaded with shape: (300, 512, 512)

MIP¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt

def plot_segmented_planes(ct_images, segmentation_images):
    # Extraer los cortes para cada plano
    axial_slice = ct_images[ct_images.shape[0] // 2, :, :]
    coronal_slice = ct_images[:, :, ct_images.shape[2] // 2]
    sagittal_slice = ct_images[:, ct_images.shape[1] // 2, :]

    # Extraer las segmentaciones para cada plano
    axial_seg = segmentation_images[ct_images.shape[0] // 2, :, :]
    coronal_seg = segmentation_images[:, :, ct_images.shape[2] // 2]
    sagittal_seg = segmentation_images[:, ct_images.shape[1] // 2, :]

    fig, axes = plt.subplots(1, 3, figsize=(18, 6))

    # Visualizar el plano Axial
    axes[0].imshow(axial_slice, cmap='gray')
    axes[0].imshow(axial_seg, cmap='viridis', alpha=0.5)
    axes[0].set_title('Axial Plane')
    axes[0].axis('off')

    # Visualizar el plano Coronal
    axes[1].imshow(coronal_slice, cmap='gray', aspect=aspect_ratio_coronal)
    axes[1].imshow(coronal_seg, cmap='viridis', alpha=0.5)
    axes[1].set_title('Coronal Plane')
    axes[1].axis('off')

    # Visualizar el plano Sagital
    axes[2].imshow(sagittal_slice, cmap='gray', aspect=aspect_ratio_sagittal)
    axes[2].imshow(sagittal_seg, cmap='viridis', alpha=0.5)
    axes[2].set_title('Sagittal Plane')
    axes[2].axis('off')

    plt.show()

plot_segmented_planes(windowed_ct_images, img_segmented)
No description has been provided for this image
In [ ]:
def calculate_mip(volume, axis):
    """Calculate the Maximum Intensity Projection (MIP) along a given axis."""
    return np.max(volume, axis=axis)

def create_mips(ct_volume, segmentation_volume):
    """Calculate MIPs for all three planes and display them with a lilac color theme, including segmentations."""
    # Calcular MIPs para las imágenes CT
    mip_axial_ct = calculate_mip(ct_volume, axis=0)
    mip_coronal_ct = calculate_mip(ct_volume, axis=1)
    mip_sagittal_ct = calculate_mip(ct_volume, axis=2)
    
    # Calcular MIPs para las segmentaciones
    mip_axial_seg = calculate_mip(segmentation_volume, axis=0)
    mip_coronal_seg = calculate_mip(segmentation_volume, axis=1)
    mip_sagittal_seg = calculate_mip(segmentation_volume, axis=2)

    return mip_axial_ct, mip_coronal_ct, mip_sagittal_ct, mip_axial_seg, mip_coronal_seg, mip_sagittal_seg

def display_mip_images(mip_axial_ct, mip_coronal_ct, mip_sagittal_ct, mip_axial_seg, mip_coronal_seg, mip_sagittal_seg):
    """Display MIP images for axial, coronal, and sagittal planes with customized color."""
    plt.figure(figsize=(18, 6))

    # Colormap personalizado que tiende al lila
    cmap = plt.get_cmap('twilight')
    cmap_seg = plt.get_cmap('viridis')  # Use a distinct colormap for segmentation overlay

    # Definir los percentiles para ajustar los valores de visualización de cada MIP
    vmin_a, vmax_a = np.percentile(mip_axial_ct, [2, 98])
    vmin_c, vmax_c = np.percentile(mip_coronal_ct, [2, 98])
    vmin_s, vmax_s = np.percentile(mip_sagittal_ct, [2, 98])

    # Axial
    plt.subplot(1, 3, 1)
    plt.imshow(mip_axial_ct, cmap=cmap, aspect='auto', vmin=vmin_a, vmax=vmax_a, origin='lower')
    plt.imshow(mip_axial_seg, cmap=cmap_seg, aspect='auto', alpha=0.7, origin='lower')
    plt.title('MIP Axial')
    plt.axis('off')

    # Coronal
    plt.subplot(1, 3, 2)
    plt.imshow(mip_coronal_ct, cmap=cmap, aspect=aspect_ratio_coronal, vmin=vmin_c, vmax=vmax_c, origin='lower')
    plt.imshow(mip_coronal_seg, cmap=cmap_seg, aspect=aspect_ratio_coronal, alpha=0.7, origin='lower')
    plt.title('MIP Coronal')
    plt.axis('off')

    # Sagittal
    plt.subplot(1, 3, 3)
    plt.imshow(mip_sagittal_ct, cmap=cmap, aspect=aspect_ratio_sagittal, vmin=vmin_s, vmax=vmax_s, origin='lower')
    plt.imshow(mip_sagittal_seg, cmap=cmap_seg, aspect=aspect_ratio_sagittal, alpha=0.7, origin='lower')
    plt.title('MIP Sagittal')
    plt.axis('off')

    plt.tight_layout()
    plt.show()
In [ ]:
mip_axial_ct, mip_coronal_ct, mip_sagittal_ct, mip_axial_seg, mip_coronal_seg, mip_sagittal_seg = create_mips(ct_images, img_segmented)
display_mip_images(mip_axial_ct, mip_coronal_ct, mip_sagittal_ct, mip_axial_seg, mip_coronal_seg, mip_sagittal_seg)
No description has been provided for this image

Animation¶

In [ ]:
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import pydicom
import scipy.ndimage

# Funciones para obtener planos MIP
def MIP_sagittal_plane(volume):
    return np.max(volume, axis=2)

def MIP_coronal_plane(volume):
    return np.max(volume, axis=1)

def rotate_on_axial_plane(volume, angle):
    """Rotar volumen en el plano axial."""
    return scipy.ndimage.rotate(volume, angle, axes=(1, 2), reshape=False)


def crop_volume(volume, start_slice, crop_slices):
    """Recortar el volumen para eliminar el ruido en los bordes."""
    return volume[:, :, start_slice:crop_slices]

def plot_3d_combined(ct_image, seg_image, pixel_len_mm, n=16, output_path='results/MIP'):
    """Plot 3D visualization of CT data and segmentation data."""
    os.makedirs(output_path, exist_ok=True)
    projections = []

    # Recortar el volumen para eliminar el ruido en los bordes
    #ct_image = crop_volume(ct_image, 50, 450)
    #seg_image = crop_volume(seg_image, 50, 450)

    # Crear y guardar las proyecciones MIP
    for idx, alpha in enumerate(np.linspace(0, 360, num=n, endpoint=False)):
        rotated_ct_img = rotate_on_axial_plane(ct_image, alpha)
        rotated_seg_img = rotate_on_axial_plane(seg_image, alpha)

        projection_ct = MIP_sagittal_plane(rotated_ct_img)
        projection_seg = MIP_sagittal_plane(rotated_seg_img)

        # Guardar las proyecciones
        projections.append((projection_ct, projection_seg))
        plt.figure(figsize=(18, 6))
        plt.imshow(projection_ct, cmap='gray', aspect=pixel_len_mm[0] / pixel_len_mm[1], alpha=0.7, origin='lower')
        plt.imshow(projection_seg, cmap=cmap, aspect=pixel_len_mm[0] / pixel_len_mm[1], alpha=0.4, origin='lower')
        plt.title(f'Projection {idx}')
        plt.axis('off')

        # Legend elements
        handles = [plt.Line2D([0], [0], color=cmap(i), lw=4) for i in range(1, 5)]
        labels = ['Liver', 'Mass', 'Portal Vein', 'Abdominal Aorta']
        plt.legend(handles, labels, loc='upper right')

        plt.savefig(f'{output_path}/Projection_{idx}.png')
        plt.close()

    # Crear la animación GIF
    fig, ax = plt.subplots()
    ani_frames = [
        [ax.imshow(proj[0], cmap='gray', aspect=pixel_len_mm[0] / pixel_len_mm[1], animated=True, alpha=0.7, origin='lower'),
         ax.imshow(proj[1], cmap=cmap,  aspect=pixel_len_mm[0] / pixel_len_mm[1], animated=True, alpha=0.4, origin='lower')]
        for proj in projections
    ]
    anim = animation.ArtistAnimation(fig, ani_frames, interval=100, blit=True)
    gif_path = os.path.join(output_path, 'Animation.gif')
    anim.save(gif_path, writer='pillow')  # Use 'pillow' instead of 'imagemagick'

    print(f"Ho he guardat a --> '{gif_path}'")
    return gif_path

# Configurar el colormap para las segmentaciones
cmap = mcolors.ListedColormap(['none', 'red', 'green', 'blue', 'yellow'])  # 'none' es para transparencia
bounds = [0, 1, 2, 3, 4, 5]
norm = mcolors.BoundaryNorm(bounds, cmap.N)

pixel_len_mm = [3.27, 0.98, 0.98] # Pixel length in mm [z, y, x]

# Crear y guardar la animación
gif_path = plot_3d_combined(ct_images, img_segmented, pixel_len_mm, n=64)

# Mostrar el GIF en el notebook
from IPython.display import Image
display(Image(filename=gif_path))
Ho he guardat a --> 'results/MIP\Animation.gif'
<IPython.core.display.Image object>
No description has been provided for this image

Image Coregistration¶

Loading and Preprocessing the Images¶

In [ ]:
from pydicom import dcmread
import pylibjpeg
from scipy.optimize import minimize
from scipy.spatial.transform import Rotation as R
import numpy as np
from scipy.ndimage import affine_transform, zoom
import matplotlib.pyplot as plt
import pydicom
import os
In [ ]:
reference_image_path = "Project/icbm_avg_152_t1_tal_nlin_symmetric_VI.dcm"
input_images_path = "Project/RM_Brain_3D-SPGR"
atlas_path = "Project/AAL3_1mm.dcm"
In [ ]:
input_images_dicom = sort_dicom_files_by_position(load_dicom_files(input_images_path))
reference_image_dicom = dcmread(reference_image_path)
atlas_image_dicom = dcmread(atlas_path)

# Convert to numpy arrays for processing
input_images = np.array([img.pixel_array for img in input_images_dicom])
reference_image = reference_image_dicom.pixel_array
atlas_image = atlas_image_dicom.pixel_array

print("Reference image shape:", reference_image.shape)
print("Input image shape:", input_images.shape)
print("Atlas image shape:", atlas_image.shape)
Reference image shape: (193, 229, 193)
Input image shape: (212, 512, 512)
Atlas image shape: (181, 217, 181)

Visualisation of the images

In [ ]:
def visualize_comparative_slices(image1, image2, slice_types=('axial', 'coronal', 'sagittal'), titles=('Reference Image', 'Input Image')):
    """
    Displays comparative slices from two images in axial, coronal, and sagittal planes.
    """
    # Indices for middle slices
    slices = {
        'axial': image1.shape[0] // 2,
        'coronal': image1.shape[1] // 2,
        'sagittal': image1.shape[2] // 2
    }

    # Setting up the figure and axes
    fig, axes = plt.subplots(nrows=2, ncols=len(slice_types), figsize=(5 * len(slice_types), 10))

    for idx, slice_type in enumerate(slice_types):
        # Reference image slice
        ax1 = axes[0, idx]
        ax1.imshow(image1[slices[slice_type], :, :] if slice_type == 'axial' else image1[:, slices[slice_type], :] if slice_type == 'coronal' else image1[:, :, slices[slice_type]], cmap='gray')
        ax1.set_title(f'{titles[0]} - {slice_type.capitalize()}')
        ax1.axis('off')

        # Input image slice
        ax2 = axes[1, idx]
        ax2.imshow(image2[slices[slice_type], :, :] if slice_type == 'axial' else image2[:, slices[slice_type], :] if slice_type == 'coronal' else image2[:, :, slices[slice_type]], cmap='gray')
        ax2.set_title(f'{titles[1]} - {slice_type.capitalize()}')
        ax2.axis('off')

    plt.tight_layout()
    plt.show()

# Usage of the function with example images
visualize_comparative_slices(reference_image, input_images)
No description has been provided for this image

Con estos resultados podemos ver que las dos imágenes tienen una orientación diferente. Para alinearlas, tendremos que tener en cuenta esta diferencia. Es por ello que vamos a aplicar un alineamiento inicial.

MATCHING DIMENSIONALITY

si feim q els dos tensors siguin 128x128, me farà que sigui idferent --> Tensor de mides diferents pero s'objecte

  • mida de cada pxel fisica (pixel sppacing --> en cada slice com te mous en x,y) slice_thinkes --> en en y,z
In [ ]:
def crop_image(img, target_shape):
    """
    Crops an image to a specified target shape by calculating the cropping indices
    directly and applying them to slice the numpy array. This version aims to use
    direct slicing for an intuitive understanding and avoids the explicit loop.
    """
    # Calculate start and end indices for each dimension to ensure the crop is centered
    start_indices = [(s - t) // 2 for s, t in zip(img.shape, target_shape)]
    end_indices = [start + t for start, t in zip(start_indices, target_shape)]

    # Prepare slice objects for each dimension
    slices = tuple(slice(start, end) for start, end in zip(start_indices, end_indices))

    # Crop the image using calculated slices
    cropped_img = img[slices]

    return cropped_img


def get_pixel_spacing_and_thickness(dicom):
    """
    Extracts PixelSpacing and SliceThickness from a DICOM file.
    Converts them from DSfloat to regular float or tuple of floats.
    """
    try:
        pixel_spacing = tuple(map(float, dicom.PixelSpacing)) if dicom.PixelSpacing else (None, None)
        slice_thickness = float(dicom.SliceThickness) if dicom.SliceThickness else None
    except AttributeError:
        print("DICOM file is missing necessary attributes.")
        pixel_spacing, slice_thickness = (None, None), None

    print("---- Header information ----")
    print(f" - PixelSpacing: {pixel_spacing}")
    print(f" - SliceThickness: {slice_thickness}")
    return pixel_spacing, slice_thickness

def match_dimension(images, shape, interpolation=cv2.INTER_LINEAR):
    """Matching dimensions with a given shape"""
        # Create an empty array to hold the resized images
    resized_stack = np.empty((images.shape[0], shape[1], shape[0]), dtype=images.dtype)

    # Resize each slice and store in the new array
    for i, img in enumerate(images):
        resized_stack[i] = cv2.resize(img, shape, interpolation=interpolation)

    return resized_stack
In [ ]:
# Stack the 2D slices into a 3D volume
input_images_volume = np.stack(input_images)

# Getting PixelSpacing and Thickness for match dimensions
print(f"\nInput Image Pixel Spacing and Slice Thickness for first image:")
pixel_spacing, slice_thickness = get_pixel_spacing_and_thickness(input_images_dicom[0])

reference_image_croped = crop_image(reference_image, atlas_image.shape) # Matching dimensions with atlas image
input_images_reshaped = match_dimension(input_images_volume,(int(input_images_volume.shape[1] * pixel_spacing[0]), int(input_images_volume.shape[2] * pixel_spacing[1])))
input_images_cropped = crop_image(input_images_reshaped, reference_image_croped.shape)

# Normalize images
reference_image_norm = reference_image_croped / np.max(reference_image_croped)
input_images_norm = input_images_cropped / np.max(input_images_cropped)

print("\nResampling images..")
print("Resampled reference image shape:", reference_image_norm.shape)
print("Resampled input image shape:", input_images_norm.shape)
print("Atlas image shape (not resampled):", atlas_image.shape)
visualize_comparative_slices(reference_image_norm, input_images_norm)
Input Image Pixel Spacing and Slice Thickness for first image:
---- Header information ----
 - PixelSpacing: (0.5078, 0.5078)
 - SliceThickness: 2.0

Resampling images..
Resampled reference image shape: (181, 217, 181)
Resampled input image shape: (181, 217, 181)
Atlas image shape (not resampled): (181, 217, 181)
No description has been provided for this image

Rigid motion & loss function¶

In [ ]:
import quaternion
import scipy.ndimage

# Compute MSE between reference image and transformed
def compute_mse(fixed, moving):
    return np.mean((fixed - moving) ** 2)

def calculate_center(image):
    """Calculate the geometric center of a 3D image."""
    return np.array(image.shape) / 2

def rigid_transform(image, params, output_shape):
    """Apply a rigid transformation using quaternions for rotation and a vector for translation."""
    if len(params) != 7:
        raise ValueError("Parameter array must have 7 elements: 4 for quaternion and 3 for translation.")

    # Unpack quaternion and translation from params
    q = np.quaternion(params[0], params[1], params[2], params[3])
    translation = np.array(params[4:7])
    
    # Normalize the quaternion to ensure it represents a rotation
    q_normalized = q / np.abs(q)
    
    # Convert the quaternion to a rotation matrix
    rotation_matrix = quaternion.as_rotation_matrix(q_normalized)
    
    # Calculate transformation offset
    center = calculate_center(image)
    rotated_center = np.dot(rotation_matrix, center)
    offset = translation - rotated_center + center

    # Perform affine transformation
    transformed_image = affine_transform(
        image,
        rotation_matrix,
        offset=offset,
        output_shape=output_shape,
        order=3,
        mode='constant',
        cval=0.0
    )
    return transformed_image

def objective_function(params, fixed_image, moving_image):
    """Calculate the MSE between the fixed image and a transformed moving image."""
    transformed_image = rigid_transform(moving_image, params, fixed_image.shape)
    return compute_mse(fixed_image, transformed_image)

def coregistration(method, options):
    """Perform image registration using optimization to minimize MSE."""
    initial_params = np.array([1, 0, 0, 0, 0, 0, 0])  # Identity quaternion and zero translation
    result = minimize(
        objective_function, 
        initial_params,
        args=(reference_image_norm, input_images_norm),
        method=method, 
        options=options
    )
    return result

def apply_transformation(image, parameters, reference_image, plot=True):
    """Apply the optimal transformation and visualize the results."""
    transformed_image = rigid_transform(image, parameters, reference_image.shape)
    if plot:
        visualize_comparative_slices(reference_image, transformed_image)
    return transformed_image

def debug_transformation(reference_image, transformed_image, title='Transformation Debug'):
    """
    Function to display side-by-side comparison to help identify issues in the transformation process.
    """
    fig, axes = plt.subplots(1, 3, figsize=(18, 6))
    axes[0].imshow(reference_image[reference_image.shape[1] // 2], cmap='gray')
    axes[0].set_title('Reference Image - Axial')
    axes[1].imshow(transformed_image[transformed_image.shape[1] // 2], cmap='gray')
    axes[1].set_title('Transformed Image - Axial')
    axes[2].imshow(np.abs(reference_image[reference_image.shape[1] // 2] - transformed_image[transformed_image.shape[1] // 2]), cmap='hot')
    axes[2].set_title('Difference Image - Axial')
    for ax in axes:
        ax.axis('off')
    plt.tight_layout()
    plt.show()


def apply_predefined_rotation(image, rotation_angles):
    """
    Apply predefined rotations to the image using quaternions.
     - image: 3D numpy array of the image to be transformed.
     - rotation_angles: Tuple of angles (in radians) to rotate around each axis (x, y, z).
    """
    # Creation the rotation matrices for rotations around the x and z axes
    rot_x = np.array([
        [1, 0, 0],
        [0, np.cos(rotation_angles[0]), -np.sin(rotation_angles[0])],
        [0, np.sin(rotation_angles[0]), np.cos(rotation_angles[0])]
    ])
    
    rot_z = np.array([
        [np.cos(rotation_angles[2]), -np.sin(rotation_angles[2]), 0],
        [np.sin(rotation_angles[2]), np.cos(rotation_angles[2]), 0],
        [0, 0, 1]
    ])

    # Convert quaternion to rotation matrix
    rotation_matrix = np.dot(rot_z, rot_x)
    
    # Calculate center and offset for transformation
    center = calculate_center(image)
    rotated_center = np.dot(rotation_matrix, center)
    offset = -rotated_center + center
    
    # Perform the transformation
    rotated_image = affine_transform(
        image, 
        rotation_matrix, 
        offset=offset, 
        output_shape=image.shape, 
        order=3, 
        mode='constant', 
        cval=0.0
    )
    return rotated_image
import matplotlib.pyplot as plt

def visualize_coronal_slices(image1, image2, title1='Reference Image - Coronal Slice', title2='Transformed Image - Coronal Slice'):
    """
    Visualize coronal plane slices of two images side by side.
    Args:
        image1 (numpy.ndarray): First image array (reference image).
        image2 (numpy.ndarray): Second image array (transformed input image).
        title1 (str): Title for the first image.
        title2 (str): Title for the second image.
    """
    # Calculate the index for the middle slice in the coronal plane
    coronal_slice_index = image1.shape[1] // 2

    # Set up the plotting environment
    fig, axes = plt.subplots(1, 2, figsize=(12, 6))

    # Plotting the coronal slice of the first image
    ax1 = axes[0]
    ax1.imshow(image1[:, coronal_slice_index, :], cmap='gray')
    ax1.set_title(title1)
    ax1.axis('off')  # Hide the axis

    # Plotting the coronal slice of the second image
    ax2 = axes[1]
    ax2.imshow(image2[:, coronal_slice_index, :], cmap='gray')
    ax2.set_title(title2)
    ax2.axis('off')  # Hide the axis

    plt.tight_layout()
    plt.show()
In [ ]:
# Apply predefined rotations to the Nelder-Mead transformed image
#rotated_nm_image = apply_predefined_rotation(input_images_norm, (np.pi, 0, np.pi))
#rotated_reference_image = apply_predefined_rotation(reference_image_norm, (np.pi, 0, np.pi))
rotated_nm_image = rigid_transform(input_images_norm, [-181,-217,0,np.pi, 0,0,0], reference_image_norm.shape)

visualize_comparative_slices(reference_image_norm, rotated_nm_image)
No description has been provided for this image

Optimization of transform parameters¶

In [ ]:
# Compute MSE before optimization
initial_mse = compute_mse(reference_image_norm, input_images_norm)
print(f"Initial MSE: {initial_mse}")

# Optimize using BFGS
options_bfgs = {'maxiter': 500, 'gtol': 1e-6}
result_bfgs = coregistration('BFGS', options_bfgs)
optimized_params_bfgs = result_bfgs.x
mse_bfgs = compute_mse(reference_image_norm, rigid_transform(input_images_norm, optimized_params_bfgs, reference_image_norm.shape))

# Optimize using Nelder-Mead
options_nm = {'maxiter': 500, 'xatol': 1e-6}
result_nm = coregistration('Nelder-Mead', options_nm)
optimized_params_nm = result_nm.x
mse_nm = compute_mse(reference_image_norm, rigid_transform(input_images_norm, optimized_params_nm, reference_image_norm.shape))

# Print results
print(f"BFGS Optimization: MSE = {mse_bfgs}, Parameters = {optimized_params_bfgs}")
print(f"Nelder-Mead Optimization: MSE = {mse_nm}, Parameters = {optimized_params_nm}")
BFGS Optimization: MSE = 0.06503201621259201, Parameters = [1. 0. 0. 0. 0. 0. 0.]
Nelder-Mead Optimization: MSE = 0.06500624308699136, Parameters = [ 1.01294041e+00 -4.51770740e-07  3.43398539e-08 -2.27624101e-07
  2.21116042e-04 -7.60686421e-05  1.16682785e-04]

Based on the results it appears that Nelder-Mead has given a better Mean Squared Error (MSE) than BFGS. Now we are going to viusailze the transformations.

In [ ]:
# Apply transformations using optimized parameters from Nelder-Mead
transformed_image_nm = rigid_transform(input_images_norm, optimized_params_nm, reference_image_norm.shape)
In [ ]:
# Apply predefined rotations to the Nelder-Mead transformed image
rotated_nm_image = apply_predefined_rotation(transformed_image_nm, (np.pi, 0, np.pi))
rotated_nm_image = rigid_transform(transformed_image_nm,)
rotated_reference_image = apply_predefined_rotation(reference_image_norm, (np.pi, 0, np.pi))

visualize_coronal_slices(rotated_reference_image, rotated_nm_image)
No description has been provided for this image
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import imageio.v2 as imageio  # Using imageio v2 explicitly
import os
from matplotlib.colors import Normalize

def create_alpha_fusion_gif(reference_img, transformed_img, filename='fusion.gif', duration=0.1):
    """
    Create an animated GIF with alpha-fusion of two 3D images (reference and transformed).
    Args:
        reference_img (numpy.ndarray): Reference image array.
        transformed_img (numpy.ndarray): Transformed image array.
        filename (str): Name of the GIF file to save.
        duration (float): Duration each frame appears in the GIF.
    """
    num_slices = reference_img.shape[1]  # Assuming slicing along the second axis (coronal)
    norm = Normalize(vmin=0, vmax=max(reference_img.max(), transformed_img.max()))

    # Create a temporary directory to store images
    temp_dir = 'temp_images'
    if not os.path.exists(temp_dir):
        os.makedirs(temp_dir)

    image_paths = []

    for i in range(num_slices):
        fig, ax = plt.subplots()
        ref_slice = reference_img[:, i, :]
        trans_slice = transformed_img[:, i, :]
        
        # Normalize and overlay the images with alpha blending
        ax.imshow(norm(ref_slice), cmap='gray', interpolation='none')
        ax.imshow(norm(trans_slice), cmap='gray', alpha=0.5, interpolation='none')  # Change alpha for blending effect
        ax.axis('off')
        
        # Save each slice as an image
        filepath = f'{temp_dir}/slice_{i}.png'
        plt.savefig(filepath, bbox_inches='tight', pad_inches=0)
        plt.close()
        
        # Append the image path to list
        image_paths.append(filepath)
    
    # Create the GIF
    images = [imageio.imread(path) for path in image_paths]
    imageio.mimsave(filename, images, duration=duration)
    
    # Clean up the temporary images
    for image_path in image_paths:
        os.remove(image_path)
    
    print(f"GIF saved as {filename}")
In [ ]:
create_alpha_fusion_gif(rotated_reference_image, rotated_nm_image)
GIF saved as fusion.gif
In [ ]:
import numpy as np
import pandas as pd

def extract_thalamus_ids(atlas_info_path):
    """ Extracts Thalamus IDs from a .txt file. """
    try:
        df = pd.read_csv(atlas_info_path, header=None, delim_whitespace=True)
        thalamus_ids = df[df[1].str.contains('Thal', na=False)][0].tolist()
    except Exception as e:
        print(f"Error reading the thalamus IDs: {e}")
        thalamus_ids = []
    return thalamus_ids

# Load thalamus IDs
thal_ids = extract_thalamus_ids('Project/AAL3_1mm.txt')
print("Thalamus IDs:", thal_ids)

# Transform the thalamus region to the input space
def transform_thalamus(atlas_img, thal_ids, parameters):
    thal_mask = np.isin(atlas_img, thal_ids)
    transformed_thal_mask = apply_transformation(thal_mask.astype(float), parameters, atlas_img.shape, plot=False)
    return transformed_thal_mask

# Ensure you have the atlas_image loaded and the optimized_params_* from your optimization routines
transformed_thalamus_bfgs = transform_thalamus(atlas_image, thal_ids, optimized_params_bfgs)
transformed_thalamus_nm = transform_thalamus(atlas_image, thal_ids, optimized_params_nm)

# Visualize the thalamus region transformation
visualize_comparative_slices(reference_image_norm, transformed_thalamus_bfgs, slice_types=('axial', 'coronal', 'sagittal'), titles=('Reference Image', 'Thalamus - BFGS Transformed'))
visualize_comparative_slices(reference_image_norm, transformed_thalamus_nm, slice_types=('axial', 'coronal', 'sagittal'), titles=('Reference Image', 'Thalamus - Nelder-Mead Transformed'))
Thalamus IDs: [81, 82, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[153], line 46
     43 optimized_params_bfgs = [1, 0, 0, 0, 0, 0, 0]  # example parameters
     45 # Transform the thalamus region to the input space
---> 46 transformed_thalamus_bfgs = apply_transformation(atlas_image, optimized_params_bfgs, atlas_image.shape)

Cell In[153], line 20, in apply_transformation(image, parameters, output_shape, plot)
     18 transformed_image = rigid_transform(image, parameters, output_shape)
     19 if plot:
---> 20     visualize_comparative_slices(output_shape, transformed_image)
     21 return transformed_image

Cell In[153], line 31, in visualize_comparative_slices(reference_image, transformed_image)
     29 """Visualize slices for comparison."""
     30 fig, ax = plt.subplots(1, 2)
---> 31 ax[0].imshow(reference_image, cmap='gray')
     32 ax[0].set_title('Reference Image')
     33 ax[1].imshow(transformed_image, cmap='gray')

File /opt/homebrew/anaconda3/envs/DL/lib/python3.10/site-packages/matplotlib/__init__.py:1465, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1462 @functools.wraps(func)
   1463 def inner(ax, *args, data=None, **kwargs):
   1464     if data is None:
-> 1465         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1467     bound = new_sig.bind(ax, *args, **kwargs)
   1468     auto_label = (bound.arguments.get(label_namer)
   1469                   or bound.kwargs.get(label_namer))

File /opt/homebrew/anaconda3/envs/DL/lib/python3.10/site-packages/matplotlib/axes/_axes.py:5751, in Axes.imshow(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, interpolation_stage, filternorm, filterrad, resample, url, **kwargs)
   5748 if aspect is not None:
   5749     self.set_aspect(aspect)
-> 5751 im.set_data(X)
   5752 im.set_alpha(alpha)
   5753 if im.get_clip_path() is None:
   5754     # image does not already have clipping set, clip to axes patch

File /opt/homebrew/anaconda3/envs/DL/lib/python3.10/site-packages/matplotlib/image.py:723, in _ImageBase.set_data(self, A)
    721 if isinstance(A, PIL.Image.Image):
    722     A = pil_to_array(A)  # Needed e.g. to apply png palette.
--> 723 self._A = self._normalize_image_array(A)
    724 self._imcache = None
    725 self.stale = True

File /opt/homebrew/anaconda3/envs/DL/lib/python3.10/site-packages/matplotlib/image.py:693, in _ImageBase._normalize_image_array(A)
    691     A = A.squeeze(-1)  # If just (M, N, 1), assume scalar and apply colormap.
    692 if not (A.ndim == 2 or A.ndim == 3 and A.shape[-1] in [3, 4]):
--> 693     raise TypeError(f"Invalid shape {A.shape} for image data")
    694 if A.ndim == 3:
    695     # If the input data has values outside the valid range (after
    696     # normalisation), we issue a warning and then clip X to the bounds
    697     # - otherwise casting wraps extreme values, hiding outliers and
    698     # making reliable interpretation impossible.
    699     high = 255 if np.issubdtype(A.dtype, np.integer) else 1

TypeError: Invalid shape (3,) for image data
No description has been provided for this image

Thalamus Region¶

In [ ]:
# Load the thalamus mask
def load_thalamus_mask(reference_image):
    thalamus_mask = np.zeros(reference_image.shape)
    thalamus_mask[30:35, 30:35, 30:35] = 1
    return thalamus_mask

# Applying the inverse transformation to the thalamus mask
def apply_inverse_transform(mask, quaternion, translation):
    center = calculate_center(mask)
    rotation = R.from_quat(quaternion)
    inverse_rotation = rotation.inv()
    
    def transform_coords(coords):
        coords_centered = coords - center[:, np.newaxis]
        rotated_coords = inverse_rotation.apply(coords_centered.T).T
        translated_coords = rotated_coords + center[:, np.newaxis] - translation[:, np.newaxis]
        return translated_coords
    
    coords = np.indices(mask.shape).reshape(3, -1)
    transformed_coords = transform_coords(coords).astype(int)
    
    # Ensure new coordinates are within bounds --> Antes me daba imágenes en negro
    valid_mask = np.all((transformed_coords >= 0) & (transformed_coords < np.array(mask.shape)[:, np.newaxis]), axis=0)
    transformed_coords = transformed_coords[:, valid_mask]
    original_coords = coords[:, valid_mask]
    
    transformed_mask = np.zeros(mask.shape, dtype=mask.dtype)
    transformed_mask[tuple(transformed_coords)] = mask[tuple(original_coords)]
    
    return transformed_mask

# Load and apply the inverse transform to the thalamus mask
thalamus_mask = load_thalamus_mask(resampled_reference_image_np)
inverse_transformed_thalamus = apply_inverse_transform(thalamus_mask, quaternion_opt, translation_opt)

Visualisation of the Thalamus Region¶

In [ ]:
# Función para visualizar la máscara del tálamo en el espacio de la imagen de entrada
def visualize_mask_in_input_space(image, mask, title='Thalamus Mask in Input Image Space'):
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    for i in range(3):
        slice_index = mask.shape[i] // 2
        if i == 0:
            slice_image = np.flipud(image[slice_index])
            slice_mask = np.flipud(mask[slice_index])
        elif i == 1:
            slice_image = np.flipud(image[:, slice_index])
            slice_mask = np.flipud(mask[:, slice_index])
        else:
            slice_image = np.fliplr(np.flipud(image[:, :, slice_index]))
            slice_mask = np.fliplr(np.flipud(mask[:, :, slice_index]))
        axes[i].imshow(slice_image, cmap='gray', alpha=0.7)
        axes[i].imshow(slice_mask, cmap='jet', alpha=0.3)
        axes[i].set_title(f"{title} - Slice {slice_index}")
        axes[i].axis('off')
    plt.tight_layout()
    plt.show()

visualize_mask_in_input_space(optimized_transformed_image, inverse_transformed_thalamus)

'''
pintar thalamus
'''
No description has been provided for this image
Out[ ]:
'\npintar thalamus\n'

Creating a gif for result visualisation

In [ ]:
# Create GIF of overlayed images
def create_thalamus_gif(reference_image, transformed_image, thalamus_mask, gif_path='coregistration_with_thalamus.gif'):
    images = []
    for i in range(reference_image.shape[0]):
        fig, ax = plt.subplots(1, 1, figsize=(5, 5))
        ax.imshow(reference_image[i], cmap='gray', alpha=0.5)
        ax.imshow(transformed_image[i], cmap='gray', alpha=0.5)
        thalamus_overlay = np.ma.masked_where(thalamus_mask[i] == 0, thalamus_mask[i])
        ax.imshow(thalamus_overlay, cmap='jet', alpha=0.5)  # Ensure the thalamus mask is clearly visible
        ax.axis('off')
        plt.title(f'Slice {i}')
        # Save the current figure as an image
        plt.savefig('temp.png', bbox_inches='tight', pad_inches=0)
        plt.close(fig)
        images.append(imageio.imread('temp.png'))
        os.remove('temp.png')
    imageio.mimsave(gif_path, images, duration=0.1)

def display_gif(gif_path):
    with open(gif_path, "rb") as file:
        display(Image(file.read()))

# Create GIF
create_thalamus_gif(resampled_reference_image_np, optimized_transformed_image, inverse_transformed_thalamus, gif_path='coregistration_with_thalamus.gif')
display_gif('coregistration_with_thalamus.gif')
/var/folders/gj/6tyw6rld1nn68w3tjknkqmkw0000gn/T/ipykernel_53640/2006559335.py:15: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  images.append(imageio.imread('temp.png'))
No description has been provided for this image

s'ha de pintar es thalamus a s'image de input a reference, ara es thalamo l'hem de conseguir a l'espai input. Invertir rotacio i despres translació.